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Abstract 


Three  different  formulations  are  presented  for  expressing  the  initial  motion  direction  of 
a  system  of  contacting  frictionless  rigid  bodies  under  gravity.  The  bodies  are  assumed  to 
have  no  initial  velocity.  The  first  formulation  expresses  the  accelerations  of  the  bodies  in 
terms  of  the  contact  forces  between  the  bodies.  The  contact  forces  are  themselves  expressed 
as  the  solution  to  a  quadratic  programming  problem.  The  second  formulation  expresses 
the  accelerations  of  the  bodies  according  to  Gauss’  “principle  of  least  constraint.”  This 
principle  is  well-known  to  apply  to  systems  with  holonomic  motion  constraints  (such  as 
joints  or  hinges);  in  this  paper,  we  show  that  the  principle  extends  to  the  nonholonomic 
constraints  that  arise  due  to  contact  between  bodies.  The  third  formulation  is  conceptually 
the  simplest;  it  says  simply  that  the  initial  acceleration  of  the  system  is  in  the  direction  that 
most  quickly  decreases  the  gravitational  potential  energy  of  the  system  without  violating 
the  contact  constraints  between  the  bodies. 


1  Introduction 


In  this  paper  we  consider  the  dynamics  of  a  collection  of  frictionless  rigid  bodies  with 
contact  constraints.  All  bodies  are  initially  motionless  and  are  acted  upon  by  an  external 
force  mg  where  m  is  a  body’s  mass  and  g  G  R3  indicates  a  gravity  field.  One  or  more  of 
the  bodies  are  assumed  to  be  fixed  in  place.  Since  the  bodies  are  initially  motionless,  the 
impending  motion  for  each  body  is  in  the  direction  of  the  initial  acceleration  of  that  body.  In 
this  paper,  we  show  the  equivalence  of  three  different  formulations  for  expressing  the  initial 
direction  of  acceleration.  The  first  two  formulations  characterize  not  only  the  direction  of 
the  initial  acceleration,  but  the  magnitude  as  well. 


The  first  formulation  for  the  acceleration  of  rigid  bodies  with  contact  constraints  has 
appeared  several  times  in  the  literature[5,  6,  8,  4,  1],  This  formulation  expresses  the 
acceleration  of  each  body  as  a  function  of  the  net  force  and  torque  acting  on  each  body. 
Since  the  external  gravity  force  mg  acting  on  each  body  is  known  a  priori,  only  the  unknown 
contact  forces  that  arise  between  bodies  at  contact  points  need  to  be  determined.  The 
contact  forces  can  be  expressed  in  terms  of  the  solution  to  a  convex  quadratic  programming 
problem.  Solving  convex  quadratic  programs  is  a  polynomial-time  problem[7]. 


The  second  formulation  is  an  application  of  Gauss’  principle  of  least  constraint  to  the 
collection  of  rigid  bodies.  Gauss’  principle  expresses  the  acceleration  of  systems  with 
holonomic  motion  constraints  as  the  solution  to  a  minimization  problem.  We  have  not 
encountered  any  application  of  Gauss’  principle  to  systems  with  nonholonomic  constraints 
in  the  literature.  We  will  show  that  Gauss’  principle  applies  to  the  nonholonomic  motion 
constraints  that  prevent  interpenetration  between  contacting  bodies.  In  the  second  formu¬ 
lation,  the  acceleration  of  the  bodies  is  expressed  as  the  solution  to  a  convex  quadratic 
programming  problem. 


The  third  formulation  is  in  some  ways  the  most  attractive  and  elegant  of  the  three.  If  we 
think  of  gravity  as  the  gradient  of  a  potential  energy  function,  we  naturally  picture  the  initial 
acceleration  of  the  system  as  a  motion  that  carries  the  system  “downhill,”  with  respect  to  the 
potential  energy  function.  We  will  show  that  the  initial  acceleration  direction  of  the  system 
is  parallel  to  the  steepest  descent  direction  down  the  potential  energy  function  that  does  not 
violate  constraints  due  to  contact.  This  appears  to  be  an  obvious  statement;  what  is  not 
so  obvious,  however,  is  that  this  statement  is  not  well-defined  until  we  describe  precisely 
how  we  measure  the  steepness  of  energy  descent  in  a  given  direction.  Such  a  definition  is 
intimately  linked  with  how  we  measure  distance  in  the  space  of  motions  for  our  system  of 
rigid  bodies. 
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2  Rigid  Body  Formulation 

2.1  Mass  Distribution 

Let  us  describe  the  mass-distribution  of  a  rigid  body  in  a  global  frame  of  reference  as  a  set 
of  mass  points,  each  with  location  p,  and  mass  m,.  The  total  mass  M  of  a  body  is 


M= 

i 


(i) 


The  vector  c  denotes  the  center  of  mass  of  the  body;  that  is,  c  satisfies 


Y  «.-(Pi  -  c)  =  o 


(2) 


for  each  body.  (We  will  denote  row  vectors,  column  vectors,  and  matrices  whose  entries 
are  all  zero  simply  by  “0”  throughout  this  paper.  The  dimension  of  0  should  be  clear  from 
the  context  in  which  it  occurs.  A  scalar  value  of  zero  is  written  simply  as  “0.”) 

If  we  let  r,  denote  the  world-space  displacement  of  the  tth  mass  point  from  the  center 
of  mass  by 

r,  =  P,  -  c  (3) 

then  the  inertia  tensor  I  of  the  body  is 


/ 

i  =  5> 

i 

\ 


nl  +  nl 


-nynx 


-rlxrt>.  ~rixriz  > 
n]  +  nj  -nynz 
-nzriy  n2x  +  r2  / 


(4) 


For  a  vector  u  G  R3,  define  U*  to  be  the  anti-symmetric  matrix 


(0  —u.  uY  \ 

uz  0  -ux  .  (5) 

-Uy  U,  0  ) 

For  any  vector  v  6  R3,  (u*)v  =  u  x  v.  Additionally,  the  relation  -u*u*  =  (uru)l  -  uur 
holds,  where  1  is  the  3  x  3  identity  matrix.  Using  this  relation,  it  is  easy  to  show  that 

I  =  £«»((r«7*v)l  -  r,r,r)  =  Y,  -”hr'r,\  (6) 

i  i 

In  section  4,  we  will  also  make  use  of  the  relations  u*v  =  —  v*u  and  u*T  =  -u*. 


2.2  Motion  Constraints 

We  will  represent  possible  motions  of  a  system  of  rigid  bodies  in  terms  of  virtual  displace¬ 
ments  of  each  body.  Let  <5p,  =  (<^c,,  1)0,)  represent  a  displacement  of  the  rth  body  in  the 
system,  with  <5c,  and  80,  vectors  in  R3.  The  vector  Ac,  denotes  a  translational  displacement 
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Figure  I:  Contact  between  bodies  A  and  B.  Motion  constraints  are  formulated  in  terms  of 
the  relative  motion  of  the  bodies  at  points  d  and  d'. 

of  the  «th  body,  while  89,  denotes  a  rotation  of  magnitude  ||60,||  of  the  body  around  its 
center  of  mass.  The  axis  of  the  rotation  is  along  the  89,  direction. 

Contact  between  bodies  generates  constraints  on  the  allowable  displacements.  Consider 
figure  1  where  bodies  A  and  B  contact.  If  body  A  undergoes  a  displacement^  =  (8c„,89„), 
then  point  d,  as  attached  to  body  A,  undergoes  a  particular  displacement  (‘id,,.  Similarly,  a 
displacement  8\th  of  body  B  causes  a  displacement  8Ab  of  point  d,  as  attached  to  body  B.  To 
prevent  interpenetration  from  occurring,  the  relative  displacement  ^d„-^d/,  cannot  have  any 
component  opposite  the  unit  normal  direction  n.  We  can  express  this  as  the  constraint 


n  •  (6d(l  -  8d„ )  >  0.  (7) 

Similarly,  we  also  need  to  prevent  interpenetration  from  occurring  at  point  d'  by  requiring 
that  n  •  (<5d^  -  <5d{,)  >  0.  If  body  B  was  fixed,  the  motion  constraint  at  d  would  simply 
be 

n  •  6d„  >  0  (8) 

and  similarly  for  d'.  To  simplify  bookkeeping,  we  do  not  count  fixed  objects  as  bodies  in 
our  system;  rather,  we  simply  note  when  regular  movable  objects  are  in  contact  with  fixed 
objects,  and  generate  the  appropriate  motion  constraint,  such  as  equation  (8). 

We  will  assume  that  the  motion  constraints  can  be  expressed  by  a  finite  number  of 
constraint  inequalities  in  the  form  of  equation  (7)  or  (8),  all  of  which  must  be  satisfied. 
(Palmer[9]  and  Baraff[2]  contain  further  discussion  on  this  issue.)  That  is,  we  consider 
systems  whose  motion  constraints  are  expressed  in  terms  of  m  contact  points  between  the 
bodies.  Let  the  ith  contact  point  of  the  system  be  a  contact  between  bodies  A  and  B  at  the 
point  d,  in  a  global  frame  of  reference.  Let  n,  denote  the  unit  surface  normal,  pointing 
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outwards  from  B  towards  A  at  d„  and  let  c„  and  ch  denote  the  positions  of  the  center  of  mass 
of  bodies  A  and  B  respectively.  If  A  undergoes  a  displacement  <5p(J  =  (be,,,  b0„)  then  d,,  as 
attached  to  A,  undergoes  the  displacement 

bcu  +  bG„  x  (d,  -  c„). 

Similarly,  for  a  displacement  ^  =  (bcb,  bdb)  of  body  B,  d,’s  displacement,  as  attached  to 
B ,  is 

bch  +  bOh  x  (d,  -  c,,). 

The  motion  constraint  at  the  ith  contact  point  is  therefore 


n,  •  (bc„  +  bO,,  x  (d,  -  c(1)  -  be i,  -  bdh  x  (d,  -  ch))  >  0.  (9) 


Since  each  constraint  is  a  linear  inequality  on  the  be  and  bO  variables,  we  can  express 
the  simultaneous  satisfaction  of  all  the  cons. faints  as  one  large  linear  system.  If  the  vector 
(*>p  denotes  the  virtual  displacements  of  the  n  bodies  by  writing 


(  W|  ) 


bp  = 


) 


then  we  express  all  m  motion  constraints  by  writing 


Ji^p  >  0 


(10) 


where  J  is  an  m  x  6 n  matrix.  The  coefficients  of  J  are  computed  according  to  equation  (9). 
Using  this  notation,  we  can  say  that  a  legal  motion  for  the  system  is  a  displacement  i^p  that 
satisfies  J6p  >  0.  Note  that  the  displacement  tfp  =  0  always  yields  a  legal  motion  (the 
null-motion). 


3  Contact  Force  Formulation 

At  each  of  the  m  contact  points,  a  contact  force  may  arise  between  the  contacting  bodies, 
to  prevent  interpenetration.  Since  we  are  dealing  with  frictionless  contacts,  we  know  that 
the  contact  forces  will  act  normal  to  the  contact  surfaces.  Thus,  at  the  rth  contact  point,  we 
consider  a  contact  force  A,n,  that  acts  on  body  A  of  the  contact,  and  a  contact  force  -A,n, 
that  acts  on  body  B  of  the  contact,  with  A,  the  unknown  scalar  magnitude  of  the  force-pair. 
Since  n,  is  directed  from  B  towards  A,  and  since  contact  forces  must  be  repulsive,  we  require 
A,  >  0  for  each  contact  point. 
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We  define  the  6/1  x  6 n  block-diagonal  generalized  mass  matrix  M  as 

/M,  1  0  -  \ 

0  1. 


M„  1  0 

V  •••  0  ij 

where  1  denotes  the  3  x  3  identity  matrix,  and  M,  and  1/  are  the  mass  and  inertia  tensor  of 
the  /th  body.  A  net  force  and  torque  of  F,  and  r,  acting  on  each  body  is  represented  as  the 
generalized  force  vector  Q  of  length  6 n,  defined  by 


The  force  exerted  on  each  body  by  gravity  is  M,g  while  the  torque  exerted  is  zero  (since  the 
gravity  field  is  uniform).  Thus,  the  generalized  gravity  force  QA,  is 


Q*  = 


(with  0g  R-  for  this  definition). 

If  v;  and  u>,  denote  the  linear  and  angular  velocity  of  the  /th  rigid  body,  the  acceleration 
a  from  a  generalized  force  Q  acting  on  the  system  is 


svr'Q. 


Let  A  e  Rm  be  the  vector  of  contact  force  magnitudes  A,.  The  generalized  force  Q,  due 
to  the  contact  force  pairs  A,n,  and  -A,n,  at  each  contact  point  is[8J 

Q,  =  JrA.  (15) 

The  acceleration  a  of  the  system  due  to  both  the  contact  forces  and  gravity  is  thus 

a  =  M-,(Qr  +  QJ  =  M  '/A  +  M-'Q,.  (16) 
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Because  our  system  is  initially  motionless,  the  constraint  JSp  >  0  on  displacements 
yields  the  constraint 

Ja  >  0  (17) 

on  the  accelerations  of  the  bodies.  Consider  the  scalar  (Ja)„  which  characterizes  the  relative 
acceleration  in  the  n,  direction  at  the  ith  contact  point.  If  this  quantity  is  positive,  contact 
is  being  broken.  Otherwise,  (Ja),  =  0  and  contact  is  not  broken.  Since  frictionless  contact 
forces  are  workless,  the  contact  force  magnitude  at  the  ith  contact  point  must  be  zero  if 
contact  is  being  broken.  We  express  this  constraint  by  writing 


A,(Ja),  =  0.  (18) 

Since  each  A,  is  nonnegative,  we  have  A  >  0.  Combining  this  with  the  constraints  Ja  >  0 
and  A,(Ja),  =  0  for  all  i,  we  can  write 

A  >  0,  Ja  >  0  and  Ar(Ja)  =  0.  (19) 

Using  equation  (16),  we  rewrite  this  as 

A  >  0,  JM'jrA  +  JM  lQ*  >  0  and  \T  (jM-'jrA  +  JM-'Q,.)  =  0.  (20) 

Equation  (20)  can  be  viewed  as  a  quadratic  program  for  the  unknown  A.  If  the  matrix 
JM~'jr  is  positive  definite  then  A  is  unique.  Otherwise  JM'J7  is  positive  semidefinite 
(since  M  and  thus  M_l  is  positive  definite),  and  a  solution  exists  for  A  although  it  may  not 
be  unique.  However,  Cottle[3]  has  shown  that  if  A|  and  A?  are  solutions  to  equation  (20) 
then 


JM-1J7Aj  =  JM'jrA2 

(21) 

which  implies 

JM-'jr(A,  -  A2)  =  0. 

(22) 

Multiplying  both  sides  of  this  equation  by  (A|  —  A2)7  yields 

(A. - 

A2)7JM-|J7(A1  -  A2)  =  (Jr(A,  -  A2))rM-|(J7(A,  -  A2))  =  0. 

(23) 

But  since  M~‘ 

is  positive  definite  it  must  be  that 

Jr(A|  -  A2)  =  0, 

(24) 

which  implies  JrA|  =  J7A2.  Thus,  even  if  the  solution  A  to  equation  (20)  is  not  unique,  the 
resulting  acceleration  of  the  system 

M-'(JrA)  +  M-'Q,, 


will  be. 

Equations  (16)  and  (20)  give  us  our  first  characterization  of  the  impending  motion  of 
the  system.  Given  any  A  that  is  a  solution  to  equation  (20),  equation  (16)  describes  the 
resulting  acceleration  a,  and  thus  the  impending  motion  direction  of  the  system. 
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4  Gauss’  Least  Constraint  Formulation 


Gauss’  principle  of  least  constraint  is  a  very  elegant  statement  of  the  acceleration  of  a  system 
with  constraints[?].  Consider  a  system  of  particles  p,  €  R\  each  with  mass  m,  and  acted 
upon  by  a  force  f,  e  R\  Let  the  scalar  quantity  Z,  called  the  “constraint”  of  the  system,  be 
defined  by 

z  =  5Z  ^(/i  -  m.P'Yif.  ~  m«P<)-  (25) 

Note  that  Z  is  nonnegative. 

Gauss’  principle  states  very  simply  that  the  accelerations  p,  of  the  particles  will  min¬ 
imize  Z.  If  all  the  particles  are  completely  unconstrained,  the  principle  is  obviously  true, 
since  the  acceleration  of  every  particle  satisfies  m,p,  =  /„  yielding  Z  =  0.  However,  if 
there  are  constraints  on  the  particles’  accelerations,  then  the  accelerations  which  minimize 
Z  subject  to  those  constraints  are  the  accelerations  that  will  actually  occur. 

Gauss’  principle  is  easily  shown  to  apply  to  systems  with  holonomic  constraints.  In 
such  a  system,  if  a  motion  direction  is  legal  the  opposite  (or  reverse)  motion  direction  is 
legal  as  well.  This  is  not  necessarily  true  in  our  system.  We  will  show  however  that  Gauss’ 
principle  can  be  applied  to  our  problem. 

Let  us  write  equation  (25)  in  the  form 

Z  =  £  £  9— (/j»  ~  mJiPji)T(fji  ~  m & ,)  (26) 

j—  i  I 

where  index  j  runs  over  the  n  bodies,  and  index  t  runs  over  the  points  of  the  yth  body.  The 
quantities  am,,,  py,  and  /y,  are  the  mass,  acceleration  and  force  acting  on  the  /th  point  of  the 
yth  body. 

From  equation  (3),  we  can  express  the  location  of  the  /th  particle  in  the  y'th  body  in  terms 
of  the y'th  body’s  center  of  mass  cy,  and  the  displacement  ry,  as 

P ji  =  9/  “h  9/1-  (22) 

The  derivative  of  a  vector  r  attached  to  a  body  with  angular  velocity  u/  is  given  by  r  =  u;  x  r, 
so  differentiating  equation  (27)  yields 

P,/  =  Vy  +  UJj  X  Tji.  (28) 

Since  vy  =  u>y  =  0,  differentiating  again  yields 

py(  =  vy  +  u>y  x  Tji.  (29) 

Using  the  “*”  notation  defined  in  section  2.1,  and  the  fact  that  u*v  =  -v‘u  we  can  rewrite 
equation  (26)  in  the  form 

"  I 

z  =  Y.Y.  ~  W  ~  x  rjfifj.  -  mjiij  -  mj^j  x  *>•) 

j=  i  I 

"  l 

=  Z  Z  nZrifp  -  "Wi  -  -  »h'*j  -  mJ<Vj  rj>)-  (3°) 

y=l  /  Lmi' 


1 


Expanding  equation  (30)  using  the  relations  d^r,,  =  —  r *tu>j  and 

=  -fJtfiUj  =  {-fJ,r*,Uj)T  =  -vjrffji  =  ufcfj, 

yields 

z  =  EE^r  (flh  ~  2mpffrj  -  2mj>flu]T]l 

j=  i  /  zmj< 


+  H^J‘  +  +  ml(u}rj>)T(vJrji)) 

7=  1  I  Lml' 

-  mjivJrJjUjj  +  \nijtf\j  +  lJmJl(r*luJ)T(r*iwJ))  . 

We  can  break  this  up  into  separate  sums: 

Z  =  EE^-t(E*)Vi>;(E«Jft 


;=1  i 
n 


rV  7=1  V  I 


j=  I 


-E*J  Ew/»rjs Uy  +  E*!  EH 


y=i 

/i 


;=i 


+  EWI  EH<r; 


7=1 


(31) 


The  following  identities  hold  for  each  body:  for  any  body  j,  the  net  force  F;  acting  on 
that  body  is 

Fj  =  E  *• 


Similarly,  the  net  torque  on  the  jth  body  is 

Tj =  E  rj‘ x  fj>  =  E  rjifj<- 

i  i 

From  equations  (1)  and  (2)  the  relations 

M}  =  E  mJ>  and  E  - 0 

hold  for  each  body,  and  using  the  linearity  of  the  operator, 

E  mi>ri  =  (Emprj>y  =  ° 

Last,  the  inertia  tensor  I,  of  each  body  is  given  by 

I;  =  ^  )  mjirjJjr 


(32) 


(33) 


(34) 


(35) 


(36) 
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Using  these  relations,  equation  (31)  simplifies  to 


Z  =  i  £  #  -  t  -  t  ufr  +  i  ±  yJUjirj  +  i  £  (37) 

7=1  i  7=1  7=1  7=1  7=1 

If  the  net  force  F;  on  each  body  is  Mjg  and  the  net  torque  is  zero,  using  the  definitions 
of  Qg,  a  and  M  from  the  previous  section  we  can  write  simply 

z = £  £  -  arQ«  +  5“rMa-  <38> 

7=1  I 

We  would  like  to  show  that  the  acceleration  a  which  solves 

minZ  subject  to  Ja  >  0  (39) 

is  the  same  as  the  acceleration  a  of  the  previous  section.  Since  Z  is  a  quadratic  function 
of  a,  and  M  is  positive  definite,  problem  (39)  is  a  convex  quadratic  programming  problem 
with  a  unique  solution  a.  The  first-order  optimality  conditions,  or  KKT  conditions!?],  for 
a  constrained  optimization  problem  are  both  necessary  and  sufficient  when  applied  to  a 
convex  quadratic  programming  problem.  The  KKT  conditions  for  a  to  be  a  solution  to 
problem  (39)  are  that  there  exists  A  €  Rm  such  that 

^  -  JrA  =  0,  A  >  0,  Ja  >  0,  and  Ar(Ja)  =  0.  (40) 

aa 

Differentiating  Z  with  respect  to  a  yields 

§  =  M»-Q„-  (41) 


The  condition  dZ/d a  —  JrA  =  0  is  thus 


Ma  -  Q„  -  JrA  =  0 


or  simply 


a  =  MlQg  +  M1JrA. 


(43) 


This  means  that  the  necessary  and  sufficient  conditions  for  a  to  solve  problem  (39)  are 
simply 

a  =  M  'Qg  +  M-lJrA,  A  >  0,  Ja  >  0,  and  Ar(Ja)  =  0.  (44) 


However,  this  is  precisely  the  same  as  the  formulation  for  a  given  by  equations  ( 1 6)  and  (20) 
in  section  3.  As  we  noted  before,  the  product  JrA  is  unique  even  though  A  may  not  be. 
Equation  (44)  is  a  direct  proof  of  this,  in  that  a  solution  a  to  equation  (44)  is  unique  because 
it  is  the  (unique)  minimizer  of  a  positive  definite  quadratic  program,  namely  problem  (39). 
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5  Steepest  Descent  Formulation 

The  formulation  of  this  section  is  the  simplest  of  the  three  to  state.  We  claim  that  a  system’s 
initial  acceleration  is  in  the  direction  which  most  quickly  the  decreases  potential  of  the 
system  without  violating  the  contact  constraints.  Although  the  physical  intuition  behind 
this  statement  is  valid,  the  translation  of  this  statement  to  a  mathematical  problem  requires 
care. 

Suppose  that  we  describe  the  geometric  state  of  the  n  bodies  in  our  system  by  a  vector 
x  e  R6n.  The  uniform  gravity  field  acting  on  our  system  of  objects  can  be  considered  the 
gradient  of  a  potential  energy  function  U(\).  We  would  like  to  describe  the  impending 
motion  of  our  system  by  saying  that  the  system  moves  in  the  direction  that  most  quickly 
decreases  the  potential  energy,  without  violating  any  of  the  contact  constraints;  that  is,  the 
system  moves  in  the  direction  of  steepest  descent  (with  respect  to  the  function  U)  that  is 
legal.  What  precisely  do  we  mean  when  we  say  “steepest  descent”? 

Given  two  displacement  directions  6p , ,  <5p2  e  R6n  we  say  that  <5p,  is  a  steeper  displace¬ 
ment  direction  if 

VU(x)  •  <5p,  <  VU(x)  •  Sp2.  (45) 

Clearly,  this  definition  makes  no  sense  unless  we  are  comparing  vectors  6p,  and  <5p2  of  equal 
length.  If  we  agree  to  define  a  motion  direction  in  our  system  as  a  displacement  <5p  of  unit 
length,  then  we  can  say  that  the  direction  of  steepest  descent  is  the  unit  displacement  <5p  that 
minimizes 

V(/(x)  •  6p. 

At  first  glance  this  appears  to  adequately  define  the  steepest  descent  direction.  In  fact, 
the  definition  is  still  incomplete.  Until  we  specify  how  we  plan  to  measure  distance,  that 
is,  what  constitutes  a  unit  vector,  we  still  cannot  say  what  the  steepest  descent  direction  is. 
For  example,  consider  figure  2  which  shows  a  potential  energy  function  U(x,y)  =  )'■  In 
figure  2a,  we  have  defined  length  using  the  two-norm  ||v]|2  of  a  vector  defined  by 

II  v||  2  =  n/^v  =  +  Vy  +  V2.  (46) 

When  we  use  this  “standard”  distance  measure,  the  set  of  unit  vectors  centered  at  a  point 
(jt,y)  in  the  plane  traces  out  a  circle  in  the  plane.  Using  this  distance  measure,  the  steepest 
descent  direction  at  (x,  y)  is 

which  is  the  direction  (0,-1)  since  VU  points  straight  upwards  (at  every  point). 

In  figure  2b  however,  we  are  using  a  different  metric  for  measuring  distance.  Given  a 
symmetric  positive  definite  matrix  A  we  define  the  “||  •  j|A-norm”  of  a  vector  v  by1 

II  v||  a  =  VvrA\.  (47) 

'The  special  case  when  A  is  the  identity  matrix  multiplied  by  some  positive  scalar  o2  yields 
llvll A  =  allvll2  f°r  v  The  matr'x  A  used  to  define  distance  in  figure  2b  is  not  a  scalar  multiple  of 
the  identity  matrix. 
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Figure  2:  (a)  The  set  of  vectors  v  defined  by  ||v||2  =  1  forms  a  circle.  The  vector  in  this  set 
for  which  the  largest  decrease  in  U  is  achieved  is  indicated  by  the  bold  vertical  vector,  (b) 
Distance  is  now  defined  by  the  ||  •  ||A-norm,  where  A  is  a  symmetric  positive  definite  matrix. 
The  set  of  vectors  v  satisfying  ||v||A  =  1  is  an  ellipse,  rather  than  a  circle.  Under  this  new 
distance  measure,  the  unit-length  vector  for  which  U  decreases  the  most  points  down  and 
to  the  left. 

Under  this  distance  metric,  the  vector  <5p  which  satisfies  ||<5pj| A  =  1  and  gives  the  steepest 
descent  (that  is,  minimizes  VU(x,y)  •  <5p)  is  not  parallel  to  VU(x,y).  Instead,  the  steepest 
descent  direction  points  down  and  to  the  left. 

Why  might  we  prefer  to  use  a  norm  other  than  the  standard  two-norm  for  measuring 
distance?  Consider  a  planar  rigid  body  with  degrees  of  freedom  x,  y.  and  9.  What  constitutes 
a  “unit  displacement”  of  the  body?  We  could  consider  the  set  of  displacements  6x,  6y  and 
66  of  the  rigid  body  such  that 

(6xf  +  (f,yf  +  (Mf=  1;  (48) 

that  is,  all  unit  two-norm  displacements.  If  we  measure  x  and  y  in  centimeters  and  9  in 
radians,  then  the  set  of  unit  displacements  is  some  set  V| .  However,  if  we  measure  x  and 
y  in  meters  and  9  in  radians,  we  get  a  completely  different  set  of  displacements  V2.  The 
sets  Vi  and  V2  are  fundamentally  different,  in  that  neither  is  a  scalar  multiple  of  the  other, 
just  as  neither  the  circle  nor  the  ellipse  in  figure  2  is  a  scalar  multiple  of  the  other.  As  a 
result,  if  we  define  the  steepest  descent  direction  based  on  a  two-norm  measure  of  distance, 
the  steepest  descent  direction  depends  on  the  units  of  measurement  chosen.  Defining  unit 
displacements  of  a  system  of  rigid  bodies  based  on  the  two-norm  is  therefore  a  completely 
arbitrary  definition. 

We  propose  that  a  more  natural  way  to  measure  the  length  of  a  displacement  is  based 


on  the  kinetic  energy  of  a  system.  If  the  vector  v  is  the  velocity  of  our  system,  that  is. 


then  the  kinetic  energy  T  of  the  system  is 

T  =  ±vTMv  (50) 

where  M  the  generalized  mass  matrix  defined  in  equation  (11).  We  will  define  distance 
in  terms  of  the  j|  •  ||M-norm;  this  allows  us  a  measurement-invariant  way  of  defining  the 
steepest  descent  direction.2 

We  would  like  to  show  that  using  the  jj  •  j|M-norm  we  can  characterize  the  initial 
acceleration  of  our  system  of  objects  in  a  very  simple  manner.  The  initial  acceleration 
direction  is  described  by  simply  saying  that  it  is  the  steepest  legal  descent  direction  down 
the  function  U,  where  steepest  is  with  respect  to  distance  measured  by  the  ||  •  ||M-norm. 
Mathematically,  we  will  say  that  if  a  system  with  configuration  x  has  nonzero  acceleration, 
then  that  acceleration  is  parallel  to  the  displacement  <5p  which  solves 

min  W(x)  •  6p  subjeci  io  j  {[*| }  •  (51) 

In  the  case  of  a  uniform  gravity  field,  —VU  =  everywhere,  so  we  can  say  that  the 
acceleration  is  parallel  to  the  solution  <5p  of 

mm  —SpTQf,  subject  to  {  ^||  !“,“}■  (52) 

(Note  that  we  do  not  bother  to  characterize  the  acceleration  direction  of  a  stable  system 
(that  is,  a  system  with  acceleration  a  =  0)  in  terms  of  a  solution  to  problem  (52).  For  such 
systems,  problem  (52)  may  fail  to  have  a  unique  solution  or  any  solution  at  all.) 

To  prove  our  claim,  we  will  show  that  whenever  the  solution  a  to  equation  (39)  is 
nonzero,  there  exists  a  positive  scalar  a  such  that 

<5p  -  aa  (53) 

2The  shape  of  the  set  of  unit  vectors  under  this  norm  does  not  depend  on  the  units  of  measurement  used. 
That  is,  consider  two  observers  with  differing  measurement  systems.  Suppose  observer  A  determines  the  set 
of  velocities  VA  for  which  the  kinetic  energy  of  the  system  is  some  particular  value  TA,  and  suppose  observer 
B  determines  the  set  of  velocities  VB  which  yield  a  kinetic  energy  of  TB.  If  the  energies  TA  and  TB  are  the 
same  (even  though  A  and  B  may  use  different  units  to  describe  them),  then  the  sets  VA  and  VB  describe  the 
same  vectors  (even  though  A’s  and  B's  coordinate  description  of  the  individual  vectors  differ,  because  A  and 
B  have  differing  measuring  systems).  Even  if  TA  and  TB  are  not  the  same,  the  sets  VA  and  VB  are  the  same 
up  to  a  scalar  multiple;  that  is,  every  vector  in  VA  corresponds  to  a  scalar  a  times  a  vector  in  VB,  where  a  is 
a  constant  depending  on  TA  and  TB.  Thus,  two  observers  can  always  select  vector  sets  with  the  same  shape 
(but  not  scale)  by  having  each  observer  select  all  vectors  that  yield  unit  kinetic  energy  (or  any  constant)  in  the 
measurement  units  chosen  by  that  observer. 
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solves  problem  (52).  Let  a  =  a*  and  A  =  A*  satisfy  equation  (44),  so  that  a*  is  the  actual 
acceleration  of  the  system,  with  a*  nonzero.  We  claim  that 


<5p  = 


1 


a* 


(54) 


is  the  solution  to  problem  (52).  Since  problem  (52)  is  a  linear  minimization  with  quadratic 
constraints,  problem  (52)’s  KKT  conditions  are  both  necessary  and  sufficient  for  a  solution 
Sp.  The  KKT  conditions  to  problem  (52)  are 

Qs  +  JrA  -  2sMSp  =  0,  ||<5p||m  =  1,  A  >  0,  J«5p  >  0,  and  Ar(J6p)  =  0  (55) 


where  s  is  an  unconstrained  scalar  and  A  €  Rm.  We  can  rewrite  these  conditions  as 


2s<5p  =  lVr'Q*  +  M_lJrA,  U^pIIm  =  1,  A  >  0,  J<5p  >  0,  and  XT(J6p)  =  0.  (56) 

To  see  that  <5p  =  a*/||a* ||m  fulfills  the  KKT  conditions  (that  is,  equation  (56)),  let 
2s  =  ||a*||M.  Then  since  a*  and  A*  satisfy  equation  (44),  choosing  A  =  A*  yields 


M-'Q,,  +  M1  JrA  =  a*  =  2s<5p 

(57) 

as  well  as 

A  >  0. 

(58) 

Since  J a*  >  0  and  A *r(Ja*)  =  0,  we  have 

Jtp  =  Jtt^  =  TnJrr  Ja*  >  o 
lla  IIm  ||a  ||m 

(59) 

and 

(60) 

Finally, 

i  1 

(61) 

We  conclude  that  the  direction  of  acceleration  is  indeed  parallel  to  the  solution  (5p  to 
problem  (52). 
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